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O ,| ABSTRACT 

, We investigate the linear theory of Kelvin-Helmholtz instability at the interface be- 

tween a partially ionized dusty outflow and the ambient material analytically. We 
model the interaction as a multifluid system in a planar geometry. The unstable modes 
are independent from the charge polarity of the dust particles. Although our results 
show a stabilizing effect for charged dust particles, the growth time scale of the grow- 
r~| | ing modes gradually becomes independent of the mass or charge of the dust particles 

O |. when the magnetic field strength increases. We show that growth time scale decreases 

' ' with increasing the magnetic field. Also, as the mass of the dust particles increases, 

the growth time scale of the unstable mode increases. 

£3 . 
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^ ! 1 INTRODUCTION 

• The Kelvin-Helmholtz (KH) instability is of much interest in 
the investigation of variety of the astrophysical phenomena 
, such as interaction between jets or outflows and the ambi- 
■^j- ■ ent medium (e.g., Watson et al. 2007; Birk & Wiechen 2002; 
O ] Rosen et al. 1999; Downes & Ray 1998; Hardee & Stone 
OO , 1997; Bodo et al. 1995; ). This instability, the simplest exam- 
' pie of shear flow instability, is a well-known phenomenon in 
fluid mechanics and astrophysics (e.g. Chandrasekhar 1961). 
• ^ • Despite the very existence of the ions, neutrals and charged 
^\ ' dust particles in jets or outflows (Markwick-Kemper, Green 
H ! & Peeters 2005; Weinberger & Armsdorfer 2004; Gueth, 
- - < Bachiller & Tafalla 2003; Shepherd 2001), the KH insta- 
bility in these systems has been studied mostly in one fluid 
approximation just for simplicity. Dusty outflows are ob- 
served in some of the starburst galaxies (e.g., Alton, Davies 
& Bianchi 1999). However, a two-fluid treatment of KH in- 
stability, taking account of the different motions of ions and 
neutrals and of the magnetic field, has been investigated by 
some authors as well (e.g., Birk et al. 2000; Watson et al. 
2004; Chhajlani & Vyas 1991; Chhajlani & Vyas 1990). 

Birk et al. (2000) studied the KH instability by taking 
into account the full incompressible dynamics of both the 
neutral and the ionized gas components with applications 
to multi-phase galactic outflow winds. In another similar 
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study, Watson et al. (2004) investigated the KH instabil- 
ity in the linear, partially ionized regime to determine its 
possible effect on entrainment in massive bipolar outflows. 
They showed that for much of the relevant parameter space, 
neutral and ions are sufficiently decoupled that the neutrals 
are unstable while ions are held in place by the magnetic 
field. Shadmehri & Downes (2007) extended this analysis 
to a layer of ions and neutrals with finite thickness. They 
showed that perturbations with wavelength comparable to 
layer's thickness are significantly affected by the thickness of 
the layer. Birk & Wiechen (2002) focused on unstable shear 
flows in partially ionized dense dusty plasma. They consid- 
ered the dust and neutral gas components so that dust and 
neutral collisions is the dominant momentum transfer mech- 
anism and dust component can interact with magnetic field 
lines, although dust charge fluctuations are negligible. They 
showed long wavelength modes can be stabilized by dust and 
neutral gas collisional momentum transfer. However, their 
analysis focused on numerical solution and multifluid nu- 
merical simulations of the problem of shear flow in weakly 
ionized, magnetized dusty plasmas, and did not concentrate 
on solving the dispersion relation for this type of flow. In 
this paper, we will obtain the dispersion relation of KH in- 
stability in a multifluid system. 

Recently, Wiechen (2006) studied the KH instability by 
doing multifluid numerical simulations in partially ionized, 
dusty plasmas for different masses and charges of dust. He 
showed a stabilizing effect for more massive dust grains with 
no dependence on the charge polarity of the dust. In this pa- 
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per, we study KH instability in a multifluid system analyti- 
cally. In particular, we apply our results to magnetized dusty 
outflows. In the next section, the basic equations and as- 
sumptions which are based on Pandey & Vladimirov (2007) 
are presented. We develop analytic estimates of the relevant 
time scales of the KH modes. In section 3, properties of the 
unstable growing modes are studied. 



2 GENERAL FORMULATION 
2.1 Basic multifluid equations 

Our basic equations and the main assumptions are similar 
to Pandey & Vladimirov (2007) . We take account of the dif- 
ferent bulk velocities and densities of the neutral, electrons, 
ions and charged dust particles on both sides of the KH 
interface. The continuity equation is 



dt 



+ V.(p J -v J -) = 0, 



(1) 



where pj and Vj are the velocity of the various plasma com- 
ponents and the neutrals, respectively. 
The momentum equations are 



= -qjUj{E H — ) - PjVjnVj, 

e,i,d 



(2) 
(3) 



E' = E + v n x B/c is the electric field in the neutral frame. 
j stands for electrons (q e = — e), ions (g; = e) and dust 
(qd = Ze), where Z is the number of charge on the grain. 
The other physical variables have their usual meanings. 
Also, the collision frequencies is (Draine et al. 1983) 



< av >i 



(4) 



where < av > jn is the rate coefficient for the momentum 
transfer by the collision of the j th particle with the neutrals: 

< av > in = 1.9 x 10" 9 cmV\ (5) 
<,7O e „=4.5xl0- 9 (^)5cmV 1 , (6) 

and for small grains, we have < av > d n~< crv >» n , but for 
grains ranging between a few Angstrom to a few microns 
(Nakano & Umebayashi 1986) 



< aV >dn 



2-8xio- B (— y- 



X(: 



a N 2 3-1 
cm s 



(7) 



' 10 _5 cm 

where a is the grain radius. In our calculation, for the ion 
mass and mean neutral mass we adopt rm = 30m p and 
m n — 2.33m p , where m p — 1.67 x 10~ 24 g is the proton 
mass. 

Defining the mass density of the bulk fluid and the bulk 
velocity as p w p n and u « v n , equations (1), (2) and (3) 
give the continuity and the momentum equations for the 
bulk fluid as 



| + V.(pu)=0, 



(8) 



/ du „ . „ J x B 
p(— + u.Vu = -VP + . 

at c 



(9) 



The next simplifying assumption is that electrons and 
ions are assumed well coupled to the magnetic field which 
implies j3 e ^> /?; 3> 1, where (3j = u> c j/vj n is the ratio of 
cyclotron Lo cj = qjB/rrijC to the collision frequencies. Based 
on this assumption and using quasi-neutrality condition, wc 
have (Pandey & Vladimirov 2007; Ciolek & Mouschovias 
1993) 



-- 1 + 6 J 

Zen d 

where 

e = [i + ^]/3i 



(10) 



(11) 



Thus, the induction equation can be written as (Pandey & 
Vladimirov 2007) 



^ = V x [(u x B) - x B] 

dt lV ; Zen d 1 



(12) 



Equations (8), (9) and (12) along with the equation 



V.B = 



(13) 



are our basic equations for Kelvin-Helmholtz instability in a 
partially ionized medium. However, we assume incompress- 
ibility for the analytical calculations. 



2.2 Linear perturbations 

For doing linear analysis, the unperturbed properties of the 
system are important. We suppose that the streaming takes 
place in the a;— direction with velocity U{z), 



U(z) 



+U for z < 
-U for z > 



where U is constant. The magnetic field is assumed parallel 
to the interface; that is, B = Be x . Finally, all unperturbed 
physical quantities are assumed constant in each medium. 

Now, we can linearize the basic equations. We perturb 
the physical variables as 

X(z,x,t) = x'(z) exp[i(ujt + k x x + k y y)]. 

Thus, 



(14) 



., / , / du' z _ 
dz 

<t>pu' x = ~k x P', 



4>pu y = -k y P' + —(k x B' y - k y B' x ), 
, , .dP' iB , dB' x 

dB' 

4>B' X = k x Bu' x + r)k x B{—j^- - ikyB' z ), 



4>B' y = k x Bu' y - r)k x B(^?- - ik x B' zh 



dB' x 
dz 

<j>B' z = k x Bu' z — r)k x Bi{k x B' y — k y B' x ), 
ik x B x -(- iky By -(- — — = 0, 



(15) 
(16) 
(17) 

(18) 

(19) 

(20) 
(21) 
(22) 
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where 



= W + k x U,T) 



c 1 + 



1 B 



4tt Zen d 4n pu) mcd ' 



(23) 



and LOmcd = (pd/p)(l/l + Q)(ZeB /m,dc) is the modified 
dust-cyclotron frequency. 

Now, we can simplify the above differential equations. 
By multiplying equation (16) by k x and equation (17) by k v 
and then adding the resulting equations, we obtain 



icj>p^ = -k 2 P' + ^f-(k x B' y - k y B' x , 



(24) 



where k 2 — k 2 + k 2 . This equation gives P' in terms of the 
other variables and the parameters. By substituting P' into 
equation (18), we obtain 

k x B , 
4tt 



p(t>( V-k 2 )u' z = ^{V-k 2 )B> z 



(25) 



where T> 



On the other hand, using equations (16) and (17) wc 



have 



k x Uy kyU x 



k x B 



{h x By kyB x ) . 
(19) and (20) •< 
{k x tiy kyU x ^) 



4irp(f) 

Also, using equation (19) and (20) we obtain 
k R' - k B' kxB ' 



(26) 



i)k x Bi 



(V-k 2 )B' z . 



(27) 



So, from equations (26) and (27) we have 



k ft' — k R' 



{rik x B/4>) 



1 - {klB 2 l^p<p) 
x(V-k 2 )B' z . 

Substituting this equation into equation (21) gives 
(rfklB 2 /^ 



i k x B i 
B, = — — u. 



1 - (klB^/^p^) 



(V-k 2 )B' z . 



(28) 



(29) 



Equations (25) and (29) are our main equation for KH anal- 
ysis. However, we can simply them into one differential equa- 
tion for u' z , 



(V - q 2 )(V - k 2 )u' z =0, 



where 

2 , 2 
q = k 



_1 , VAk x . 2[ , 



2 -If, 



(30) 
(31) 



4:Tvpri 2 v <f> ' Ly VAk i 
and va = B/y/A-Kp is Alfven speed. 

2.3 Dispersion relation 

The general solution of equation (30) is a linear combination 
of exp(±fcz) and exp(±g«). Thus, we can write the solutions 
of equation (30) appropriate to the two regions as 

,,,f A ie +kz + A 2 e +qiz Iotz<Q 
Uz(Z) ~ \ A 3 e- kz + A 4 e- q * z for z > 

where A\, A2, A3 and A4 are the constants of integration yet 
to be obtained from the boundary conditions appropriate to 



our system. Note that qi and q-z are given by equation (31), 
and are assumed to have a positive real part so as to render 
the perturbations bounded at infinity. We have 



q f = k 2 {l-( 



UJ m cd,2,M 2 {x+l) 2 



kU 



x + 1 

<2/ iV2 



-] 2 }, 



(32) 
(33) 



2 ,2 M ,U mcd , 2,M (l-l) Z -l ]2l 
q2=k{1 -^ [ —1 1 } 

and M = U/va and x = uj/k x U. 

The boundary conditions to be satisfied at the interface 
z = are the standard conditions which have been used 
by many authors (Hunter & Whitaker 1989; Roychoudhury 
& Lovelace 1986; Mehta & Bhatia 1988; Chhajlani & Vyas 
1990; Chhajlani & Vyas 1991). We are using the continuity 
of normal and tangential components of the magnetic field, 
the vertical displacement of the interface and the total pres- 
sure as the boundary conditions. After doing mathematical 
manipulations, we obtain 

Ai + A 2 A 3 + A 4 



x + 1 



[1-M 2 (x + 1) 2 ]A 2 = 



x + 1 
x-1 



[1 -M 2 (x- 1) 2 ]A 4 , 



x + 1 



(kAi + qi A 2 ) = -kA 3 - q 2 A 4 , 



fcAi + giX 2 (a; + l) 2 A2 



x+1 



x-1 



x[-kA a -q 2 M 2 {x-l) 2 A 4 }. 



(34) 



(35) 



(36) 



(37) 



For a non-trivial solution of the above algebraic equations, 
the determinant of the coefficients should vanish. We get the 
following equation 



3a qi +q 2 = 



(l-/3)(« 2 -l) 



rfc, 



(38) 



(1 + /?)[1-M 2 (a; + 1) 2 ] 

where a = (x + l)/(x - 1) and (3 = [1 - M 2 (x - 1) 2 ]/[1 - 
M 2 (x + l) 2 ]. Equation (38) is the general dispersion rela- 
tion. However, it would be extremely unwieldy to solve it 
its present form. After substituting the values of qi and q 2 
and doing long algebraic mathematical manipulations, wc 
can transform equation (38) into a polynomial form, i.e. 

P S 2 8 + P 7 Z 7 + P 6 Z & + P 5 Z 5 + P 4 Z 4 + P 3 Z 3 + P 2 Z 2 



+Piz + P = 0, 



(39) 



where x = +^/z/M and the coefficients Pi ... P% are pre- 
sented in the Appendix. We note that not all roots of the 
equation (39) are acceptable. All roots are inserted into our 
fundamental dispersion relation (38) to determine which are 
the valid solutions. We accept only those solutions which 
have negative imaginary part and give us positive values for 
the real parts of q\ and q 2 . 



3 ANALYSIS 

Now, we carry out a parameter study of the roots of equa- 
tion (39) as functions of the input parameters. Although 
this dispersion relation is applicable to a wide range of as- 
trophysical systems, we restrict our analysis to dusty out- 
flows with small grains in order to illustrate the effect of 
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Figure 1. Growth time scales of the instability vs. wavelength 
of the perturbations for different level of the magnetic strength 
and various sizes of dust particles. In these plots, it is assumed 
\Z\ = 1. Curves of top plot are labeled by the size of the dust 
particles in Angstrom. But since for B > 50/iG, the curves are not 
very sensitive to the size of the dust particles or their electrical 
charge, curves of the bottom plot are labeled by the magnetic 
strength only. 



charged dust particles on the growth rate of the unstable 
modes. We take the molecular hydrogen density to be 10 3 
cm 3 and the ratio of ionized to neutral mass density is 1CP 6 . 
Also, the ratio of dust to neutral density is assumed to be 
0.01. The bulk density of a dust particle is assumed 1.25 
gr cm~ 3 . Also, the relative velocity 2U is assumed to be 
20 Km s _1 . We note that the electrical charge of the dust 
particles appears only through the modified dust-cyclotron 
frequency u> m cd- But in our final dispersion relation, this pa- 
rameter / = (u>mcd/kU)~ 2 appears and so, the polarity of 
the electrical charge of the dust particle is not important. 

Since we fix the densities and the velocity for the flow, 
when the magnetic field strength increases the Mach number 
decreases. Another important input parameter is modified 
dust-cyclotron frequency u> mc d which is calculated based on 
the other input parameters. Figure 1 shows growth rates of 
the fastest unstable growing mode versus wavelength of the 
perturbations for different level of the magnetic strength and 
various sizes of dust particles. Here, the electrical charge of 
dust particles is assumed to be \Z\ = +1. Top plot of Figure 
1 shows the dispersion relation for weak magnetic field, but 
the bottom plot is for magnetic field with stronger strength. 
We noticed as the strength of the magnetic field increases, 
the time-scale of the growing modes becomes independent of 



the variations of mass of the dust particles. In top plot of this 
figure, solid curves are corresponding to magnetic strength 
B — 10 fiG, but dashed lines are for B = 20 ^iG. Each curve 
is labeled by the size of the dust particles. When the size 
(i.e., mass) of the dust particles increases, the time-scale of 
the growing modes increases too. In other words, mass of the 
dust particles has a stabilizing effect on the growing modes. 

In bottom plot of Figure 1, the dispersion relation is 
shown for stronger magnetic field strengths, i.e. B = 60, 
80 and 100 /iG. However, the plots are insensitive to the 
mass of small dust particles. We see as the magnetic field 
increases, the growth rate of the unstable modes at a partic- 
ular wavelength of the perturbation increases. It means that 
the magnetic field has a destabilizing effect on the unstable 
modes. 

Figure 2 shows the dispersion relation for the magnetic 
strength B = 10 yitG and dust particles with radius 100 
Angstrom, but changing the electrical charge of the grains. 
As the charge of the grain increases, the growing time-scale 
of the unstable modes increases. 

We can understand the latter result as follows. From 
the form of our dispersion relation it is clear that the only 
variable which changes as we vary the grain radius, a, is / 
through the modified dust cyclotron frequency, io mc d- Using 
relations (4), (7), the definition of to mcd and recalling that 
the mass of a dust grain is proportional to a 3 we can see 
that as a increases co mc d decreases. Physically this means 
that larger grains are less well coupled to the magnetic field, 
since the dust grains can only respond to variations of the 
magnetic field with frequency less than uj mc d- As Figure 1 
shows, the magnetic field destabilizes the flow. Hence it is no 
surprise that as the dust component of the flow becomes less 
coupled to the field the flow becomes more stable. A simi- 
lar argument explains the results in Fig 2. According to the 
definition of the modified dust cyclotron frequency, we have 
Umcd oc \Z\/Q. Since equation (11) shows that oc \Z\ 2 , wc 
can conclude that uj mc d oc Thus, for increasing val- 

ues of \Z\ the dust cyclotron frequency uj m cd decreases which 
implies less magnetic coupling and greater growth times. In- 
terestingly, numerical simulations (e.g. Wiechen 2006) and 
laboratory experiments (e.g. Luo, Angelo & Merlino 2001) 
have also shown that the Kelvin-Helmholtz mode is signifi- 
cantly stabilized with increasing mass of the dust. 



4 DISCUSSION 

We studied the KH instability in a multifluid system ana- 
lytically. Our results show that a higher mass of the dust 
particles has a stabilizing effect on the growing KH modes. 
But when the strength of the magnetic field increases, the 
time-scale of the growing modes gradually are becoming in- 
dependent of dust size or mass. Also, the dynamics of unsta- 
ble modes are independent from the charge polarity of the 
dust. Considering the wavelength of the perturbations and 
the time scale of the growing modes, we may conclude that 
some of structures which are seen in dusty outflows (e.g., 
Weinberger & Armsdorfer 2004; Gueth, Bachiller & Tafalla 
2003; Shepherd 2001) are produced by the KH instability. 
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Figure 2. Growth time scales of the instability vs. wavelength of 
the perturbations for different electrical charge of the dust parti- 
cles. 
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THE COEFFICIENTS OF THE DISPERSION 
EQUATION 

The coefficients of equation (39) arc 



M 



-8M 4 



4M 8 + 28M 4 + 4M 2 /, 



Ps 
Pr 

P 6 = -8M' 

P 5 = -20MV - 56M" + 48M D + 8M*f + 8M S , 
Pa 

+36M 8 + 8M 10 + 6M 12 - 120M 6 
P 3 



-4M 6 f + 70M 4 - 20M 4 / + 4/ 2 + 40M 2 / 



8 f + 32 M 10 + 160 M e + 8 M 12 — 40 M 2 f 



+32 M 4 f - 144 M 8 - 56 M 4 + 40 M & f - 16 M 8 f, 
P 2 = 8 M 14 + 8 M 2 / 2 + 40 M 8 / - 120 M 6 - 4 M 10 f 
-8 M 4 / 2 + 4 f + 36 M 12 + 196 M s - 4 M 16 + 28 M 4 
-40 M 4 / - 144 M 10 - 16 M 6 / + 20 Af 2 /, 
Pr = -120 M s + 8 M 2 / 2 + 8 fM 12 - 20 M 10 / 
-120 M 12 -4M 2 / + 48M 6 + 8M 4 / 2 + 32M 8 / 
+24 M 4 f — 8 M 4 + 48 M 14 + 160 M 10 - 8 M 16 
—40 M 6 /, 

Po = -8 M 1S + M 20 - 20 fM 12 + 4 M 14 f - 8 M 6 
+28 M 8 + 4 M 4 / 2 + 28 M 16 + 4 f 2 M s - 40 M 8 / 



+40 M 10 / - 4 M 4 f + 70 M 12 - 8 f 2 M 6 
-56 M w + M 4 + 20 M 6 /. 
where / = (uj mcd /kU)~ 2 . 



56 M 
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